Concrete Compressive Strength¶
Concrete is the most important material in civil engineering. The concrete compressive strength is a highly nonlinear function of age and ingredients. These ingredients include cement, blast furnace slag, fly ash, water, superplasticizer, coarse aggregate, and fine aggregate.
The actual concrete compressive strength (MPa) for a given mixture under a specific age (days) was determined from laboratory. Data is in raw form (not scaled).
Variable Information:¶
Given is the variable name, variable type, the measurement unit and a brief description. The concrete compressive strength is the regression problem. The order of this listing corresponds to the order of numerals along the rows of the database.
Name -- Data Type -- Measurement -- Description
- Cement (component 1) -- quantitative -- kg in a m3 mixture -- Input Variable
- Blast Furnace Slag (component 2) -- quantitative -- kg in a m3 mixture -- Input Variable
- Fly Ash (component 3) -- quantitative -- kg in a m3 mixture -- Input Variable
- Water (component 4) -- quantitative -- kg in a m3 mixture -- Input Variable
- Superplasticizer (component 5) -- quantitative -- kg in a m3 mixture -- Input Variable
- Coarse Aggregate (component 6) -- quantitative -- kg in a m3 mixture -- Input Variable
- Fine Aggregate (component 7) -- quantitative -- kg in a m3 mixture -- Input Variable
- Age -- quantitative -- Day (1~365) -- Input Variable (Concrete hardens with time and strength increases. Usually, concrete is tested after 28 days.)
- Concrete compressive strength -- quantitative -- MPa -- Output Variable
Let's frame the problem:¶
- It is clearly a typical supervised learning task since you are given labeled training examples (each instance comes with the expected output, i.e., the Concrete compressive strength ).
- It is also a typical regression task, since you are asked to predict a value. More specifically, this is a multivariate regression problem since the system will use multiple features to make a predictionSelect a Performance Measure¶
Next step is to select a performance measure. A typical performance measure for regression problems is the Root Mean Square Error (RMSE). It gives an idea of how much error the system typically makes in its predictions, with a higher weight for large errors.
Import necessary modules¶
In [1]:
#Import all the necessary modules
import pandas as pd
import numpy as np
import os
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams["figure.figsize"] = (10,8)
import warnings
warnings.filterwarnings('ignore')
#!jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10
from scipy.stats import zscore
from sklearn.decomposition import PCA
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_nodeinteractivity = 'all'
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor, AdaBoostRegressor
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.svm import SVR
import xgboost
Load the Data¶
In [2]:
df = pd.read_csv("concrete_1.csv")
Take a Quick Look at the Data Structure¶
In [3]:
df.head()
Out[3]:
In [4]:
#The info() method is useful to get a quick description of the data, in particular the total number of rows, and each attribute’s type and number of non-null values
df.info()
In [5]:
#The describe() method shows a summary of the numerical attributes
df.describe().T
Out[5]:
In [6]:
def basic_details(df):
b = pd.DataFrame()
b['Missing value'] = df.isnull().sum()
b['N unique value'] = df.nunique()
b['dtype'] = df.dtypes
return b
basic_details(df)
Out[6]:
Observations¶
- There are 9 columns and 1030 rows
- There are no missing values in this dataset. However, slag, ash and superplastic columns have minimum values as zero's, need to check if they indicate missing values
- All the columns are numeric columns, there are no categorical columns in this dataset
- The std row shows the standard deviation, which measures how dispersed the values are
- The 25%, 50%, and 75% rows show the corresponding percentiles: a percentile indicates the value below which a given percentage of observations in a group of observations falls. For example, 25% of the observation have water component lower than 164, while 50% have water compoenent lower than 185 and 75% have water compoenent lower than 192. These are often called the 25th percentile (or 1st quartile), the median, and the 75th percentile (or 3rd quartile).
In [7]:
target = 'strength'
X = df.loc[:, df.columns!=target]
y = df.loc[:, df.columns==target]
Split the data into train and test¶
While spliting the data, we may want to ensure that the test set is representative of the various categories of strength in the whole dataset. Since the strength is a continuous numerical attribute, we will have to trick Python into interpreting your continuous numerical y variable as a categorical variable instead. By creating bins, and passing y variable into an ndarray containing those bins and the corresponding values.
In [8]:
# Create the bins. My `y` variable has 1030 observations, and I want 50 bins.
bins = np.linspace(0, df.shape[0], 50)
# Save your Y values in a new ndarray,
# broken down by the bins created above.
y_binned = np.digitize(y, bins)
# Pass y_binned to the stratify argument,
# and sklearn will handle the rest
#Test train split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y_binned,random_state=2)
Feature Selection Using Random Forest¶
Discover and Visualize the Data to Gain Insights¶
We have taken a quick glance at the data to get a general understanding of the kind of data. Now the goal is to go a little bit more in depth. First, we will make sure we have put the test set aside and we are only exploring the training set. Also, if the training set is very large, we may want to sample an exploration set, to make manipulations easy and fast. In our case, the set is quite small so we can just work directly on the full set.
In [9]:
df_temp = pd.concat([X_train,y_train],axis=1)
Lets look at the distribution plot (Univariate Analysis)¶
In [10]:
df_temp.hist(bins=30, figsize=(20,15))
plt.show()
In [11]:
import itertools
cols = [i for i in df_temp.columns if i not in 'strength']
length = len(cols)
cs = ["b","r","g","c","m","k","lime","c"]
fig = plt.figure(figsize=(13,25))
for i,j,k in itertools.zip_longest(cols,range(length),cs):
plt.subplot(5,2,j+1)
ax = sns.distplot(df_temp[i],color=k,rug=True)
ax.set_facecolor("w")
quartile_1,quartile_3 = np.percentile(df[i],[25,75])
plt.axvline(df_temp[i].median(),linestyle="dashed",label="median",color="k")
plt.axvline(df_temp[i].mean(),linestyle="dashed",label="mean",color="b")
plt.axvline(np.percentile(df_temp[i],25),linestyle="dashed",label="q1",color="r")
plt.axvline(np.percentile(df_temp[i],75),linestyle="dashed",label="q3",color="g")
plt.legend(loc="best")
plt.title(i,color="navy")
plt.xlabel("")
In [12]:
ax = sns.distplot(df_temp["strength"],color=k,rug=True)
ax.set_facecolor("w")
plt.axvline(df_temp["strength"].median(),linestyle="dashed",label="median",color="k")
plt.axvline(df_temp["strength"].mean(),linestyle="dashed",label="mean",color="b")
plt.axvline(np.percentile(df_temp["strength"],25),linestyle="dashed",label="q1",color="r")
plt.axvline(np.percentile(df_temp["strength"],75),linestyle="dashed",label="q3",color="g")
plt.legend(loc="best")
plt.title(i,color="navy")
plt.xlabel("")
Out[12]:
Observations¶
- Distribution of strength,cement variable is close to normal.
- However, we can see some long tails toward right for slag, water ,age and superplastic components.
Plot strength against age¶
In [13]:
sns.lineplot("age","strength",data=df)
Out[13]:
We see that with age the strength of cement increases upto a certain level, and then after 100 days there is a sharp decrease in the strength and also after 250 days there seems to be decrease in the strength, it could be due to more moisture absorption from the atmosphere.
In [14]:
sns.set_style('darkgrid')
g = sns.FacetGrid(df_temp,hue="age",palette='coolwarm',size=6,aspect=2)
g = g.map(plt.hist,'strength',bins=20,alpha=0.7)
In [15]:
df_temp["age"].value_counts().sort_values()
Out[15]:
Box Plots to visualize the outliers¶
In [16]:
sns.set(context="paper", font="monospace")
# Create a figure instance
fig = plt.figure(1, figsize=(18, 12))
# Create an axes instance
ax = fig.add_subplot(111)
g = sns.boxplot(data=df_temp, ax=ax, color="blue")
g.set_xticklabels(df_temp.columns,rotation=90)
# Add transparency to colors
for patch in g.artists:
r, g, b, a = patch.get_facecolor()
patch.set_facecolor((r, g, b, .3))
Observations¶
- There are outliers within slag, age, water component. The outliers present in numeric feature can be dealt by following ways based on the domain knowledge
- delete the outliers (there could be loss of data due to this)
- impute ouliers with mean/median
- impute with lower and upper bound values.
For now, we would impute the ouliers with 1% and 95% values
In [17]:
def outlier(df,columns):
for i in columns:
quartile_1,quartile_3 = np.percentile(df[i],[25,75])
quartile_f,quartile_l = np.percentile(df[i],[1,95])
IQR = quartile_3-quartile_1
lower_bound = quartile_1 - (1.5*IQR)
upper_bound = quartile_3 + (1.5*IQR)
print(i,lower_bound,upper_bound,quartile_f,quartile_l)
df[i].loc[df[i] < lower_bound] = quartile_f
df[i].loc[df[i] > upper_bound] = quartile_l
outlier(df_temp, ['age','fineagg','water'])
In [18]:
def correlation_matrix(df):
corrmat = df.corr()
top_corr_features = corrmat.index
plt.figure(figsize=(10,8))
#plot heat map
g=sns.heatmap(df[top_corr_features].corr(),annot=True,cmap="RdYlGn")
correlation_matrix(df_temp)
In [19]:
# Lets check for highly correlated variables
cor= df_temp.corr()
cor.loc[:,:] = np.tril(cor,k=-1)
cor=cor.stack()
cor[(cor > 0.55) | (cor< -0.55)]
Out[19]:
Pair Plot¶
In [20]:
sns.pairplot(df_temp,diag_kind='kde')
Out[20]:
Observations¶
- We can see density plots along the diagonal, there are multiple gaussians in every parameter
- cement parameter seems to have around 4 gaussians (3 are clearly visible and 1 is slightly visible)
- slag parameter seems to have 2 gaussians
- ash parameter seems to have 2 gaussians, first gaussian contains mostly zero's
- similarly, we can see 4 gaussians for water, superplastic, coarseagg, fineagg
- age seems to have around 5 gaussians
- Cement has a good correlation with strength, the more cement stronger the concrete.
- With age, we can see that the strength of cement decreases, it could be due to more moisture absorption from the atmosphere.
- Also, we can see that higher the water, the strength of concrete decreases, its good to have lesser water
- superplasticizers has negative correlation with water, superplasticizers always help decrease water - hence an increase in strength
- Following Features are listed based on the order of importance (ascending order) to predict strength
- cement
- age
- superplastic
- water
- fineagg
- coarseagg
- slag
- ash
We will analyse these features further, to identify key features to predict stength
Blast Furnace Slag¶
- Blast furnace slag powder as partial replacement of cement, from pair panel we can see that there is a negative correlation between cement and slag components i.e as the proportion of slag component increases the proportion of cement is low.
- Also, having zero values in slag column seems to be acceptable as its a partial replacement of cement and in cases where the values are zero we can consider that slag component was not being used.
Fly Ash¶
- Fly Ash, is also used in replacement of cement, from pair panel agaian we can see that there is negative correlation between cement and Fly Ash component
- Even here having zero values in Fly Ash, column seems to be acceptable as its a partial replacement of cement and in cases where the values are zero we can consider that Fly Ash component was not being used.
- Fly Ash has negative correlation with water
In [21]:
## check Blast Furnace Slag, Ash and Cement proportion
df_temp[(df_temp['slag']!=0) & (df_temp['ash']!=0) ].head()
Out[21]:
Superplasticizer¶
- Plasticizers are also often used when Fly ash is added to concrete to improve strength.
- Superplasticizers always help decrease batching water thus porosity - hence an increase in strength
- Having zero values in Superplasticizers also column seems to be acceptable as its a admixture used to improve the strength of concrete and not a required component.
In [22]:
df_temp[(df_temp['ash']!=0) & (df_temp['superplastic']==0)]
Out[22]:
In [23]:
## Function to train and evaluate models
def train_evaluate_model(X_train,y_train,X_test,y_test):
classifiers = [
Ridge(alpha=0.5),
Lasso(alpha=0.1),
SVR(),
GradientBoostingRegressor(random_state=2),
AdaBoostRegressor(random_state=2),
RandomForestRegressor(random_state=2),
LinearRegression(),
xgboost.XGBRegressor(random_state=2),
DecisionTreeRegressor(random_state=2)
]
# Logging for Visual Comparison
log_cols=["Classifier", "RMSE Train","RMSE Test", "R2 Train", "R2 Test"]
log = pd.DataFrame(columns=log_cols)
for clf in classifiers:
clf.fit(X_train, y_train)
name = clf.__class__.__name__
y_pred_train = clf.predict(X_train)
mse_train = np.sqrt(mean_squared_error(y_train, y_pred_train))
r2_train = r2_score(y_train, y_pred_train)
y_pred_test = clf.predict(X_test)
mse_test = np.sqrt(mean_squared_error(y_test, y_pred_test))
r2_test = r2_score(y_test, y_pred_test)
log_entry = pd.DataFrame([[name,mse_train,mse_test,r2_train,r2_test]], columns=log_cols)
log = log.append(log_entry)
log.set_index(["Classifier"],inplace=True)
return log.sort_values(by=['RMSE Test'])
In [24]:
initial_model_asis = train_evaluate_model(X_train,y_train,X_test,y_test)
initial_model_asis
Out[24]:
Observation¶
- We see that Linear models are not giving a good performance, mainly because from the pair panel we can see that there is no linear relationship between the independent variables and the target (strength).
Lets try creating polynomial Features using degree=2 and check if that boosts the performance¶
In [25]:
from sklearn.preprocessing import PolynomialFeatures
polynomial_features = PolynomialFeatures(degree=2)
x_poly_train = polynomial_features.fit_transform(X_train)
x_poly_test = polynomial_features.fit_transform(X_test)
model_poly = train_evaluate_model(x_poly_train,y_train,x_poly_test,y_test)
model_poly
Out[25]:
Observation:¶
- With polynomial of degree 2 we see that there is a boost in performance for Linear models
- SVR is performing very poorly. So we will exclude SVR this model from further steps
Feature Selection Using Random Forest¶
In [26]:
fs_model = RandomForestRegressor(n_estimators=100,oob_score=True,random_state=42)
fs_model.fit(X_train, y_train)
fs_model.feature_importances_
pd.DataFrame(fs_model.feature_importances_, columns = ["Imp"], index = X_train.columns).sort_values(['Imp'],ascending=False).plot(kind='barh',figsize=[7,6])
Out[26]:
Feature Selection Using DecisionTreeRegressor¶
In [27]:
fs_model_dt = DecisionTreeRegressor(random_state=42)
fs_model_dt.fit(X_train, y_train)
fs_model_dt.feature_importances_
pd.DataFrame(fs_model_dt.feature_importances_, columns = ["Imp"], index = X_train.columns).sort_values(['Imp'],ascending=False).plot(kind='barh',figsize=[7,6])
Out[27]:
Inference¶
- With both RandomForestRegressor and DecisionTreeRegressor, we see that cement, age, superplastic, slag and water are the key features
Feature Engineering¶
A key factor in concrete strength is the water-cement ratio, so lets create a new feature names "water_cement_ratio"
In [28]:
df_temp["water_cement_ratio"] = df_temp["water"]/df_temp["cement"]
Next, we would create "coarse_fine_agg_ratio" using Coarse and Fine aggregate components
In [29]:
df_temp["coarse_fine_agg_ratio"] = df_temp["coarseagg"]/df_temp["fineagg"]
With a bit of ground work around the domain , I also found that water_binder_ratio is another factor which could help in predicting the cement strength. So, lets create a new feature - 'water_binder_ratio'.
In [30]:
df_temp["water_binder_ratio"] = df_temp["water"]/(df_temp["cement"] + df_temp["ash"] + df_temp["slag"])
In [31]:
df_temp.head()
Out[31]:
In [32]:
sns.pairplot(df_temp,diag_kind='kde')
Out[32]:
In [33]:
correlation_matrix(df_temp)
We see a lot of correlation with the composite features and the base features using which they are derived. So. lets drop cement, water , coarseagg , fineagg , ash, slag columns as we have derived new features out of these¶
In [34]:
# #Drop a column
df_temp.drop('cement', axis=1, inplace=True)
df_temp.drop('water', axis=1, inplace=True)
In [35]:
df_temp.drop('coarseagg', axis=1, inplace=True)
df_temp.drop('fineagg', axis=1, inplace=True)
In [36]:
df_temp.drop('ash', axis=1, inplace=True)
df_temp.drop('slag', axis=1, inplace=True)
In [37]:
df_temp.info()
Lets get pre- process the test data to create the new features¶
In [38]:
df_temp_test = X_test.copy()
df_temp_test["water_cement_ratio"] = df_temp_test["water"]/df_temp_test["cement"]
df_temp_test ["coarse_fine_agg_ratio"] = df_temp_test ["coarseagg"]/df_temp_test ["fineagg"]
df_temp_test["water_binder_ratio"] = df_temp_test["water"]/(df_temp_test["cement"] + df_temp_test["ash"] + df_temp_test["slag"])
# #Drop a column
df_temp_test.drop('cement', axis=1, inplace=True)
df_temp_test.drop('water', axis=1, inplace=True)
df_temp_test.drop('coarseagg', axis=1, inplace=True)
df_temp_test.drop('fineagg', axis=1, inplace=True)
df_temp_test.drop('ash', axis=1, inplace=True)
df_temp_test.drop('slag', axis=1, inplace=True)
Build a Model with the composite features created using feature engineering¶
In [39]:
target = 'strength'
X_train = df_temp.loc[:, df_temp.columns!=target]
y_train = df_temp.loc[:, df_temp.columns==target]
X_test = df_temp_test.loc[:, df_temp_test.columns!=target]
fe_model1 = train_evaluate_model(X_train,y_train,X_test,y_test)
fe_model1
Out[39]:
In [40]:
fs_model = RandomForestRegressor(n_estimators=100,oob_score=True,random_state=42)
fs_model.fit(X_train, y_train)
fs_model.feature_importances_
pd.DataFrame(fs_model.feature_importances_, columns = ["Imp"], index = X_train.columns).sort_values(['Imp'],ascending=False).plot(kind='barh',figsize=[7,6])
Out[40]:
In [41]:
fs_model_dt = DecisionTreeRegressor(random_state=42)
fs_model_dt.fit(X_train, y_train)
fs_model_dt.feature_importances_
pd.DataFrame(fs_model_dt.feature_importances_, columns = ["Imp"], index = X_train.columns).sort_values(['Imp'],ascending=False).plot(kind='barh',figsize=[7,6])
Out[41]:
Lets check if quadratic features can help improve the model¶
In [42]:
from sklearn.preprocessing import PolynomialFeatures
polynomial_features = PolynomialFeatures(degree=2)
x_poly_train = polynomial_features.fit_transform(X_train)
x_poly_test = polynomial_features.fit_transform(X_test)
fe_model_poly = train_evaluate_model(x_poly_train,y_train,x_poly_test,y_test)
fe_model_poly
Out[42]:
Observation¶
- Linear models are not giving good performance, even with polynomial features of degree 2 for the compposite features. So, we will not consider Linear models.
- Decision Tree, RandomForest and Adaboost seems to be highly overfitting the data.
Let us explore the data for any clusters as there appears to be multiple gaussians in the independent features¶
In [43]:
from sklearn.cluster import KMeans
#check the optimal k value
ks = range(1, 9)
inertias = []
for k in ks:
model = KMeans(n_clusters=k)
model.fit(df_temp.drop('strength',axis=1))
inertias.append(model.inertia_)
plt.figure(figsize=(8,5))
plt.style.use('bmh')
plt.plot(ks, inertias, '-o')
plt.xlabel('Number of clusters, k')
plt.ylabel('Inertia')
plt.xticks(ks)
plt.show()
In [44]:
from scipy.stats import zscore
# df_std = df_temp.drop(['strength'],axis=1).apply(zscore)
df_std = df_temp.apply(zscore)
In [45]:
kmeans = KMeans(n_clusters=3,random_state=4)
kmeans.fit(df_std)
Out[45]:
In [46]:
# Check the number of data in each cluster
labels = kmeans.labels_
counts = np.bincount(labels[labels>=0])
print(counts)
In [47]:
df_std.columns
Out[47]:
In [48]:
# Distribution looks fine.
# let us check the centers in each group
centroids = kmeans.cluster_centers_
centroid_df = pd.DataFrame(centroids, columns = list(df_temp))
centroid_df.transpose()
Out[48]:
In [49]:
predictions = kmeans.predict(df_std)
predictions
df_std["group"] = predictions
df_std['group'] = df_std['group'].astype('category')
df_std.dtypes
Out[49]:
In [50]:
# Visualize the centers
df_std.boxplot(by = 'group', layout=(3,4), figsize=(15, 10))
Out[50]:
In [51]:
# Addressing outliers at group level
def replace(group):
median, std = group.median(), group.std() #Get the median and the standard deviation of every group
outliers = (group - median).abs() > 2*std # Subtract median from every member of each group. Take absolute values > 2std
group[outliers] = group.median()
return group
data_corrected = (df_std.groupby('group').transform(replace))
concat_data = data_corrected.join(pd.DataFrame(df_std['group']))
In [52]:
concat_data.boxplot(by = 'group', layout=(2,4), figsize=(15, 10))
Out[52]:
Note: When we remove outliers and replace with median or mean, the distribution shape changes, the standard deviation becomes tighter creating new outliers. The new outliers would be much closer to the centre than original outliers so we accept them without modifying them
No distinct clusters seems to be visible, lets explore further.¶
In [53]:
df_group0=df_std[(df_std['group']==0)]
df_group1=df_std[(df_std['group']==1)]
df_group2=df_std[(df_std['group']==2)]
In [54]:
sns.pairplot(df_group0,diag_kind='kde')
Out[54]:
In [55]:
sns.pairplot(df_group1,diag_kind='kde')
Out[55]:
In [56]:
sns.pairplot(df_group2,diag_kind='kde')
Out[56]:
In [57]:
# ash vs strength
var = 'water_binder_ratio'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=concat_data,hue='group')
#plot.set(ylim = (-3,3))
In [58]:
# age vs strength
var = 'age'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=concat_data,hue='group')
In [59]:
# superplastic vs strength
var = 'superplastic'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=concat_data,hue='group')
In [60]:
# water_cement_ratio vs strength
var = 'water_cement_ratio'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=concat_data,hue='group')
In [61]:
# coarse_fine_agg_ratio vs strength
var = 'coarse_fine_agg_ratio'
with sns.axes_style("white"):
plot = sns.lmplot(var,'strength',data=concat_data,hue='group')
Observation¶
There don't seem to be a clear seperation of clusters among features. It seems to be unlikely to get better results by splitting the data into clusters.
Lets tune GradientBoost and Xgboost models further to get a better performance¶
Set the range of hyperparameters to be tuned with RandonSearch for GradientBoost. Hyperparameters are the parameters that are initialized before training a model because they cannot be learned from the algorithm. They control the behavior of the training algorithm and have a high impact on the performance of a model.
In [62]:
from sklearn.model_selection import RandomizedSearchCV
import pprint
#Hyperparameters are the parameters that are initialized before training a model
#because they cannot be learned from the algorithm. They control the behavior of the training algorithm
#and have a high impact on the performance of a model
learning_rate = [0.01, 0.5, 0.1,1]
# Number of trees in GradientBoost
n_estimators = [int(x) for x in np.linspace(start = 100 , stop = 1000, num = 4)] # returns evenly spaced 10 numbers
# Number of features to consider at every split
max_features = ['auto', 'sqrt']
# Maximum number of levels in tree
max_depth = [int(x) for x in np.linspace(3, 10, num = 4)] # returns evenly spaced numbers can be changed to any
max_depth.append(None)
#the fraction of observations to be used in individual tree
subsample = [0.8,0.9,1]
# Minimum number of samples required to split a node
min_samples_split = [0.5, 1.0, 2, 5, 8, 10]
# Minimum number of samples required at each leaf node
min_samples_leaf = [1, 2, 4]
# Create the random grid
rs_param_grid = {'n_estimators': n_estimators,
'max_features': max_features,
'max_depth': max_depth,
'min_samples_split': min_samples_split,
'min_samples_leaf': min_samples_leaf,
'learning_rate':learning_rate,
'subsample':subsample}
pprint.pprint(rs_param_grid)
Tune GBM using randomisedsearchcv¶
In [63]:
gbr_clf = GradientBoostingRegressor(random_state=2)
rf_random = RandomizedSearchCV(estimator=gbr_clf, param_distributions=rs_param_grid,
n_iter = 10, scoring='neg_mean_squared_error',
cv = 5, verbose=2, random_state=42, n_jobs=-1,
return_train_score=True)
# Fit the random search model
rf_random.fit(X_train, y_train)
Out[63]:
In [64]:
rf_random.best_params_
Out[64]:
Evaluate the score of GBM by RandomisedSearchCV¶
In [65]:
best_random = rf_random.best_estimator_
best_random.score(X_test , y_test)
Out[65]:
Set the hyperparameters in the range closer to the one output by RandomizedSearchCV, to tune using GrinSearchCV
In [66]:
from sklearn.model_selection import GridSearchCV
param_grid = {
'subsample': [0.8,0.85,0.9,0.95,1],
'n_estimators': [800,900,1000,1200],
'min_samples_split': [2,4,6],
'min_samples_leaf': [2],
'max_features': ['auto','sqrt'],
'max_depth': [3,5],
'learning_rate': [0.15, 0.1,0.05,0.01]
}
Tune the GBM model using GridSearchCV¶
In [67]:
grid_search = GridSearchCV(estimator = gbr_clf, param_grid = param_grid,
cv = 5, n_jobs = -1, verbose = 2, return_train_score=True,scoring='neg_mean_squared_error')
In [68]:
# Fit the grid search to the data
grid_search.fit(X_train, y_train);
In [69]:
grid_search.best_params_
Out[69]:
Evaluate the score for GridSerachCV¶
In [70]:
best_grid = grid_search.best_estimator_
best_score = best_grid.score(X_test, y_test)
best_score
Out[70]:
In [71]:
best_grid
Out[71]:
Save the model to a file, for later use.
In [72]:
import pickle
filename = 'concrete_model.sav'
pickle.dump(best_grid,open(filename,'wb'))
Let's try to tune Xgboost parameter and see if we get better performance than GBM¶
In [73]:
xgb_model = xgboost.XGBRegressor(learning_rate=0.01,
colsample_bytree = 1,
subsample = 0.8,
objective='reg:squarederror',
n_estimators=1000,
max_depth=3,
gamma=1)
In [74]:
eval_set = [(X_train, y_train), (X_test, y_test)]
%time xgb_model.fit(X_train, y_train, eval_metric='rmse', eval_set=eval_set, verbose=True)
evals_result = xgb_model.evals_result()
RandomizedSearchCV to tune XGB parameters¶
In [75]:
learning_rate = [0.05, 0.10, 0.15, 0.20, 0.25, 0.30 ]
# Number of trees in GradientBoost
n_estimators = [int(x) for x in np.linspace(start = 100 , stop = 1000, num = 4)] # returns evenly spaced 10 numbers
# Maximum number of levels in tree
max_depth = [int(x) for x in np.linspace(3, 10, num = 6)] # returns evenly spaced numbers can be changed to any
#the fraction of observations to be used in individual tree
subsample = [0.8,0.9,1]
colsample_bytree = [0.8,0.9,1]
gamma = [0.0, 0.1, 0.2 , 0.3, 0.4, 1, 5, 10]
objective=['reg:squarederror']
#
#reg_alpha = [0.3,0.4,0.5]
# Create the random grid
rs_param_grid_xgb = {'n_estimators': n_estimators,
'max_depth': max_depth,
'subsample': subsample,
'colsample_bytree': colsample_bytree,
'gamma': gamma,
'learning_rate':learning_rate,
'objective':objective
}
pprint.pprint(rs_param_grid_xgb)
In [76]:
rf_random_xgb = RandomizedSearchCV(estimator=xgb_model, param_distributions=rs_param_grid_xgb,
n_iter = 10, scoring='neg_mean_squared_error',
cv = 5, verbose=2, random_state=42, n_jobs=-1,
return_train_score=True)
# Fit the random search model
rf_random_xgb.fit(X_train, y_train)
Out[76]:
In [77]:
rf_random_xgb.best_params_
Out[77]:
In [78]:
best_random_xgb = rf_random_xgb.best_estimator_
best_random_xgb.score(X_test , y_test)
Out[78]:
GridSearch using XGB¶
In [79]:
param_grid_xgb = {
'colsample_bytree': [0.8, 0.9, 1],
'subsample': [0.8,0.9,1],
'n_estimators': [600,800,1000],
'max_depth': [3, 5, 7, 10],
'gamma': [0, 1],
'learning_rate': [0.01, 0.05, 0.1],
#'reg_alpha': [0.3, 0.4, 0.5],
'objective': ['reg:squarederror']
}
In [80]:
grid_search_xgb = GridSearchCV(estimator = xgb_model, param_grid = param_grid_xgb,
cv = 5, n_jobs = -1, verbose = 2, return_train_score=True,scoring='neg_mean_squared_error')
In [81]:
# Fit the grid search to the data
grid_search_xgb.fit(X_train, y_train);
In [82]:
grid_search_xgb.best_params_
Out[82]:
In [83]:
best_grid_xgb = grid_search_xgb.best_estimator_
best_score_xgb = best_grid_xgb.score(X_test, y_test)
best_score_xgb
Out[83]:
GradientBoost and XGBoost are giving the same performance, i.e around 92% of accuracy score. Lets use Bootstraping technique to give the confidence interval of the model¶
In [85]:
from sklearn.utils import resample
import pickle
filename = 'concrete_model.sav'
load_lr_model=pickle.load(open(filename,'rb'))
df_train = X_train.join(y_train)
df_test = X_test.join(y_test)
df_row_merged = pd.concat([df_train, df_test], ignore_index=True)
values = df_row_merged.values
#values = df_test.values
# configure bootstrap
n_iterations = 1000
n_size = int(len(X_train) * 1)
# run bootstrap
stats = list()
for i in range(n_iterations):
# prepare train and test sets
train = resample(values, n_samples=n_size)
test = np.array([x for x in values if x.tolist() not in train.tolist()])
# fit model
model = load_lr_model
model.fit(train[:,:-1], train[:,-1])
y_test_array = test[:,-1]
# evaluate model
predictions = model.predict(test[:,:-1])
score = model.score(test[:, :-1] , y_test_array)
stats.append(score)
# plot scores
plt.hist(stats)
plt.show()
# confidence intervals
alpha = 0.95
p = ((1.0-alpha)/2.0) * 100
lower = max(0.0, np.percentile(stats, p))
p = (alpha+((1.0-alpha)/2.0)) * 100
upper = min(1.0, np.percentile(stats, p))
print('%.1f confidence interval %.1f%% and %.1f%%' % (alpha*100, lower*100, upper*100))
Conclusion¶
- We were given a dataset to predict the concrete compressive strength, given the independent features - cement, water, fly ash, slag, superplastic , age.
- We found that this is a regression problem and choose RMSE as the metric to evaluate the performance of the model
- Then we performed univariate, bivariate analysis to understand the data and found that not all the variables were highly correlated with the target.
- There were few outliers within slag, age, water component, we imputed the ouliers with 1% and 95% values
Then, created composite features -
- water_cement_ratio
- coarse_fine_agg_ratio
- water_binder_ratio
We observered a lot of correlation with the composite features and the base features using which they are derived. So. lets drop cement, water , coarseagg , fineagg , ash, slag columns as we have derived new features out of these
We also explored the data for any clusters as there appeared to be multiple gaussians in the independent features. We split the data into 3 clusters using KMean clustering. However, no distinct clusters seemed to be visible with the avaiable features.
Then we trained the model using various algorithms and evalutaed the performance of the model
- Linear models were not giving good performance, even with polynomial features of degree 2 for the compposite features.So, we did not consider Linear models.
- Decision Tree, RandomForest and Adaboost seemed to be highly overfitting the data.
- We considered GradientBoost and XGBoost for further tuning
- Linear models were not giving good performance, even with polynomial features of degree 2 for the compposite features.So, we did not consider Linear models.
We tuned the hyperparameters for GradientBoost and XGBoost using RandomizedSearch and GridSearch, found that both were giving similar performance.
At 95% confidence interval our model performs between - 88.2 % and 92.6%¶
In [ ]: